/*
* gpu_util.c
*
*  Created on: Jun 21, 2013
*      Author: tifuera
*/
#include <stdio.h>
#include <stdlib.h>
#include "gpu_util.h"

void initCublasHandle(cublasHandle_t* handle)
{
	cublasStatus_t stat;
	stat = cublasCreate(handle);
	if (stat != CUBLAS_STATUS_SUCCESS)
	{
		printf("CUBLAS INITIALIZATION FAILED");
	}

}

void copyBandMatrixFromDevice(gpu_symm_band_matrix gpu_matrix, double* h_matrix,
	cublasHandle_t handle)
{
	int num_tiles;
	int bs = gpu_matrix.block_size;
	int order = gpu_matrix.order;
	int hb = gpu_matrix.half_bandwith;
	int i, j;
	int cur_row;
	int cur_bs;
	int index;
	int cur_len;
	int cur_col;
	cublasStatus_t status;
	double* temp_tile;

	num_tiles = (order + bs - 1) / bs;
	temp_tile = (double*) malloc( (bs + hb) * bs *sizeof(double) );

	for (i = 0, cur_row = 0; i < num_tiles; i++, cur_row += bs)
	{
		cur_bs = bs;
		if (cur_row + cur_bs > order)
		{
			cur_bs = order - cur_row;
		}

		status = cublasGetMatrix(cur_bs, gpu_matrix.tile_len[i], sizeof(double), gpu_matrix.gpu_matrix_tiles[i], cur_bs, temp_tile, bs);
		if(status != CUBLAS_STATUS_SUCCESS)
		{
			printf("COPY MATRIX FROM DEVICE FAILED %d", status);
		}


		for (j = 0; j < cur_bs; j++)
		{
			index = get_index(cur_row + j, cur_row, hb, order);
			/*cublasGetVector(j + 1, sizeof(double),
				&gpu_matrix.gpu_matrix_tiles[i][j * cur_bs], 1,
				&h_matrix[index], 1);*/
			for(int it = 0; it < j +1; it++)
			{
				h_matrix[index + it] = temp_tile[j*bs + it];

			}
		}

		for (j = cur_bs; j < min(hb +1, gpu_matrix.tile_len[i]); j++)
		{
			index = get_index(cur_row + j, cur_row, hb, order);
			/*cublasGetVector(cur_bs, sizeof(double),
				&gpu_matrix.gpu_matrix_tiles[i][j * cur_bs], 1,
				&h_matrix[index], 1);*/
			for(int it = 0; it < cur_bs; it++)
			{
				h_matrix[index + it] = temp_tile[j*bs + it];

			}
		}

		cur_len = cur_bs - 1;
		int cur_col = cur_row + 1;
		for (j = hb + 1; j < gpu_matrix.tile_len[i]; j++)
		{
			index = get_index(cur_row + j, cur_col, hb, order);

			/*cublasGetVector(cur_len, sizeof(double),
				&gpu_matrix.gpu_matrix_tiles[i][j * cur_bs + (j - hb)], 1,
				&h_matrix[index], 1);*/

			for(int it = 0; it < cur_len; it++)
			{
				h_matrix[index + it] = temp_tile[j*bs + it + (j -hb)];
			}
			
			cur_col++;
			cur_len--;
		}
	}
	free(temp_tile);
}

void copyBandMatrixToDevice(double *h_matrix, gpu_symm_band_matrix * gpu_matrix,
	cublasHandle_t handle)
{
	int num_tiles;
	int bs = gpu_matrix->block_size;
	int order = gpu_matrix->order;
	int hb = gpu_matrix->half_bandwith;
	int i, j, l;
	int cur_row = 0;
	int tile_end;
	int index;
	int tile_width;
	int cur_bs;
	cublasStatus_t status;
	double* block_row;
	double* temp_tile;

	num_tiles = (order + bs - 1) / bs;
	gpu_matrix->gpu_matrix_tiles = (double **) malloc(
		num_tiles * sizeof(double *));
	gpu_matrix->tile_len = (int *) malloc(num_tiles * sizeof(int));
	block_row = (double *) malloc(bs*sizeof(double));

	temp_tile = (double*) malloc( (bs + hb) * bs *sizeof(double) );

	for (i = 0; i < num_tiles; i++)
	{
		tile_end = cur_row + bs - 1 + hb;
		if (tile_end >= order)
		{
			tile_end = order - 1;
		}
		cur_bs = bs;
		if (cur_row + bs > order)
		{
			cur_bs = order - cur_row;
		}

		gpu_matrix->tile_len[i] = tile_end - cur_row + 1;
		checkCudaErrors(
			cudaMalloc(&(gpu_matrix->gpu_matrix_tiles[i]), cur_bs * gpu_matrix->tile_len[i] * sizeof(double)));

		for (l = 0; l < cur_bs; l++)
		{
			index = get_index(cur_row + l, cur_row, hb, order);

			for(int ti = 0; ti< l+1; ti++)
			{
				temp_tile[l*bs + ti] = h_matrix[index + ti];
			}
			for(int ti = l+1; ti < cur_bs; ti++)
			{
				temp_tile[l*bs + ti] = 0;
			}
		}


		//copy full rows
		for (l = cur_bs; l < min(gpu_matrix->tile_len[i], hb +1); l++)
		{
			index = get_index(cur_row + l, cur_row, hb, order);
			for(int ti =0; ti < cur_bs; ti++)
			{

				temp_tile[l*bs + ti] = h_matrix[index + ti];
			}
		}

		//copy lower triangle
		for (l = hb + 1, j = cur_bs - 1; l < gpu_matrix->tile_len[i]; l++, j--)
		{

			index = get_index(cur_row + l, cur_row + (bs - j), hb, order);
			int q= cur_bs - j;
			int t = 0;
			
			for(int ti =0; ti < bs; ti++)
			{

				temp_tile[l*bs + ti] = 0;
			}
			for(int ti = cur_bs - j; ti < cur_bs; ti++, t++)
			{

				temp_tile[l*bs + ti] = h_matrix[index + t];
			}
		}
		status = cublasSetMatrix(cur_bs, gpu_matrix->tile_len[i], sizeof(double), temp_tile, bs, gpu_matrix->gpu_matrix_tiles[i], cur_bs);
		if(status != CUBLAS_STATUS_SUCCESS)
		{
			printf("COPY MATRIX TO DEVICE FAILED %d", status);
		}

		cur_row += bs;
	}
	printf("\nMatrix is copied to GPU\n");
	free(temp_tile);
}
